home *** CD-ROM | disk | FTP | other *** search
/ Sprite 1984 - 1993 / Sprite 1984 - 1993.iso / src / lib / m / trig.h < prev   
C/C++ Source or Header  |  1988-07-11  |  9KB  |  224 lines
  1. /*
  2.  * Copyright (c) 1987 Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms are permitted
  6.  * provided that this notice is preserved and that due credit is given
  7.  * to the University of California at Berkeley. The name of the University
  8.  * may not be used to endorse or promote products derived from this
  9.  * software without specific prior written permission. This software
  10.  * is provided ``as is'' without express or implied warranty.
  11.  *
  12.  * All recipients should regard themselves as participants in an ongoing
  13.  * research project and hence should feel obligated to report their
  14.  * experiences (good or bad) with these elementary function codes, using
  15.  * the sendbug(8) program, to the authors.
  16.  *
  17.  *    @(#)trig.h    5.2 (Berkeley) 4/29/88
  18.  */
  19.  
  20. #if defined(vax)||defined(tahoe)
  21. #ifdef vax
  22. #define _0x(A,B)    0x/**/A/**/B
  23. #else    /* vax */
  24. #define _0x(A,B)    0x/**/B/**/A
  25. #endif    /* vax */
  26. /*thresh =  2.6117239648121182150E-1    , Hex  2^ -1   *  .85B8636B026EA0 */
  27. /*PIo4   =  7.8539816339744830676E-1    , Hex  2^  0   *  .C90FDAA22168C2 */
  28. /*PIo2   =  1.5707963267948966135E0     , Hex  2^  1   *  .C90FDAA22168C2 */
  29. /*PI3o4  =  2.3561944901923449203E0     , Hex  2^  2   *  .96CBE3F9990E92 */
  30. /*PI     =  3.1415926535897932270E0     , Hex  2^  2   *  .C90FDAA22168C2 */
  31. /*PI2    =  6.2831853071795864540E0     ; Hex  2^  3   *  .C90FDAA22168C2 */
  32. static long threshx[]    = { _0x(b863,3f85), _0x(6ea0,6b02)};
  33. static long PIo4x[]    = { _0x(0fda,4049), _0x(68c2,a221)};
  34. static long PIo2x[]    = { _0x(0fda,40c9), _0x(68c2,a221)};
  35. static long PI3o4x[]    = { _0x(cbe3,4116), _0x(0e92,f999)};
  36. static long PIx[]    = { _0x(0fda,4149), _0x(68c2,a221)};
  37. static long PI2x[]    = { _0x(0fda,41c9), _0x(68c2,a221)};
  38. #define thresh    (*(double*)threshx)
  39. #define PIo4    (*(double*)PIo4x)
  40. #define PIo2    (*(double*)PIo2x)
  41. #define PI3o4    (*(double*)PI3o4x)
  42. #define PI    (*(double*)PIx)
  43. #define PI2    (*(double*)PI2x)
  44. #else   /* defined(vax)||defined(tahoe) */
  45. static double
  46. thresh    =  2.6117239648121182150E-1    , /*Hex  2^ -2   *  1.0B70C6D604DD4 */
  47. PIo4    =  7.8539816339744827900E-1    , /*Hex  2^ -1   *  1.921FB54442D18 */
  48. PIo2    =  1.5707963267948965580E0     , /*Hex  2^  0   *  1.921FB54442D18 */
  49. PI3o4    =  2.3561944901923448370E0     , /*Hex  2^  1   *  1.2D97C7F3321D2 */
  50. PI    =  3.1415926535897931160E0     , /*Hex  2^  1   *  1.921FB54442D18 */
  51. PI2    =  6.2831853071795862320E0     ; /*Hex  2^  2   *  1.921FB54442D18 */
  52. #ifdef national
  53. static long fmaxx[]    = { 0xffffffff, 0x7fefffff};
  54. #define   fmax    (*(double*)fmaxx)
  55. #endif    /* national */
  56. #endif    /* defined(vax)||defined(tahoe) */
  57. static double
  58.     zero = 0,
  59.     one = 1,
  60.     negone = -1,
  61.     half = 1.0/2.0, 
  62.     small = 1E-10,    /* 1+small**2 == 1; better values for small:
  63.              *        small    = 1.5E-9 for VAX D
  64.              *            = 1.2E-8 for IEEE Double
  65.              *            = 2.8E-10 for IEEE Extended
  66.              */
  67.     big = 1E20;    /* big := 1/(small**2) */
  68.  
  69. /* sin__S(x*x) ... re-implemented as a macro
  70.  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
  71.  * STATIC KERNEL FUNCTION OF SIN(X), COS(X), AND TAN(X) 
  72.  * CODED IN C BY K.C. NG, 1/21/85; 
  73.  * REVISED BY K.C. NG on 8/13/85.
  74.  *
  75.  *        sin(x*k) - x
  76.  * RETURN  --------------- on [-PI/4,PI/4] , where k=pi/PI, PI is the rounded
  77.  *                x    
  78.  * value of pi in machine precision:
  79.  *
  80.  *    Decimal:
  81.  *        pi = 3.141592653589793 23846264338327 ..... 
  82.  *    53 bits   PI = 3.141592653589793 115997963 ..... ,
  83.  *    56 bits   PI = 3.141592653589793 227020265 ..... ,  
  84.  *
  85.  *    Hexadecimal:
  86.  *        pi = 3.243F6A8885A308D313198A2E....
  87.  *    53 bits   PI = 3.243F6A8885A30  =  2 * 1.921FB54442D18
  88.  *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2    
  89.  *
  90.  * Method:
  91.  *    1. Let z=x*x. Create a polynomial approximation to 
  92.  *        (sin(k*x)-x)/x  =  z*(S0 + S1*z^1 + ... + S5*z^5).
  93.  *    Then
  94.  *      sin__S(x*x) = z*(S0 + S1*z^1 + ... + S5*z^5)
  95.  *
  96.  *    The coefficient S's are obtained by a special Remez algorithm.
  97.  *
  98.  * Accuracy:
  99.  *    In the absence of rounding error, the approximation has absolute error 
  100.  *    less than 2**(-61.11) for VAX D FORMAT, 2**(-57.45) for IEEE DOUBLE. 
  101.  *
  102.  * Constants:
  103.  * The hexadecimal values are the intended ones for the following constants.
  104.  * The decimal values may be used, provided that the compiler will convert
  105.  * from decimal to binary accurately enough to produce the hexadecimal values
  106.  * shown.
  107.  *
  108.  */
  109.  
  110. #if defined(vax)||defined(tahoe)
  111. /*S0     = -1.6666666666666646660E-1    , Hex  2^ -2   * -.AAAAAAAAAAAA71 */
  112. /*S1     =  8.3333333333297230413E-3    , Hex  2^ -6   *  .8888888888477F */
  113. /*S2     = -1.9841269838362403710E-4    , Hex  2^-12   * -.D00D00CF8A1057 */
  114. /*S3     =  2.7557318019967078930E-6    , Hex  2^-18   *  .B8EF1CA326BEDC */
  115. /*S4     = -2.5051841873876551398E-8    , Hex  2^-25   * -.D73195374CE1D3 */
  116. /*S5     =  1.6028995389845827653E-10   , Hex  2^-32   *  .B03D9C6D26CCCC */
  117. /*S6     = -6.2723499671769283121E-13   ; Hex  2^-40   * -.B08D0B7561EA82 */
  118. static long S0x[]    = { _0x(aaaa,bf2a), _0x(aa71,aaaa)};
  119. static long S1x[]    = { _0x(8888,3d08), _0x(477f,8888)};
  120. static long S2x[]    = { _0x(0d00,ba50), _0x(1057,cf8a)};
  121. static long S3x[]    = { _0x(ef1c,3738), _0x(bedc,a326)};
  122. static long S4x[]    = { _0x(3195,b3d7), _0x(e1d3,374c)};
  123. static long S5x[]    = { _0x(3d9c,3030), _0x(cccc,6d26)};
  124. static long S6x[]    = { _0x(8d0b,ac30), _0x(ea82,7561)};
  125. #define S0    (*(double*)S0x)
  126. #define S1    (*(double*)S1x)
  127. #define S2    (*(double*)S2x)
  128. #define S3    (*(double*)S3x)
  129. #define S4    (*(double*)S4x)
  130. #define S5    (*(double*)S5x)
  131. #define S6    (*(double*)S6x)
  132. #else    /* IEEE double  */
  133. static double
  134. S0     = -1.6666666666666463126E-1    , /*Hex  2^ -3   * -1.555555555550C */
  135. S1     =  8.3333333332992771264E-3    , /*Hex  2^ -7   *  1.111111110C461 */
  136. S2     = -1.9841269816180999116E-4    , /*Hex  2^-13   * -1.A01A019746345 */
  137. S3     =  2.7557309793219876880E-6    , /*Hex  2^-19   *  1.71DE3209CDCD9 */
  138. S4     = -2.5050225177523807003E-8    , /*Hex  2^-26   * -1.AE5C0E319A4EF */
  139. S5     =  1.5868926979889205164E-10   ; /*Hex  2^-33   *  1.5CF61DF672B13 */
  140. #endif
  141.  
  142. #if defined(vax)||defined(tahoe)
  143. #define sin__S(z)    (z*(S0+z*(S1+z*(S2+z*(S3+z*(S4+z*(S5+z*S6)))))))
  144. #else     /* defined(vax)||defined(tahoe) */
  145. #define sin__S(z)    (z*(S0+z*(S1+z*(S2+z*(S3+z*(S4+z*S5))))))
  146. #endif     /* defined(vax)||defined(tahoe) */
  147.  
  148. /* cos__C(x*x) ... re-implemented as a macro
  149.  * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS)
  150.  * STATIC KERNEL FUNCTION OF SIN(X), COS(X), AND TAN(X) 
  151.  * CODED IN C BY K.C. NG, 1/21/85; 
  152.  * REVISED BY K.C. NG on 8/13/85.
  153.  *
  154.  *                   x*x    
  155.  * RETURN   cos(k*x) - 1 + ----- on [-PI/4,PI/4],  where k = pi/PI,
  156.  *                   2    
  157.  * PI is the rounded value of pi in machine precision :
  158.  *
  159.  *    Decimal:
  160.  *        pi = 3.141592653589793 23846264338327 ..... 
  161.  *    53 bits   PI = 3.141592653589793 115997963 ..... ,
  162.  *    56 bits   PI = 3.141592653589793 227020265 ..... ,  
  163.  *
  164.  *    Hexadecimal:
  165.  *        pi = 3.243F6A8885A308D313198A2E....
  166.  *    53 bits   PI = 3.243F6A8885A30  =  2 * 1.921FB54442D18
  167.  *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2    
  168.  *
  169.  *
  170.  * Method:
  171.  *    1. Let z=x*x. Create a polynomial approximation to 
  172.  *        cos(k*x)-1+z/2  =  z*z*(C0 + C1*z^1 + ... + C5*z^5)
  173.  *    then
  174.  *      cos__C(z) =  z*z*(C0 + C1*z^1 + ... + C5*z^5)
  175.  *
  176.  *    The coefficient C's are obtained by a special Remez algorithm.
  177.  *
  178.  * Accuracy:
  179.  *    In the absence of rounding error, the approximation has absolute error 
  180.  *    less than 2**(-64) for VAX D FORMAT, 2**(-58.3) for IEEE DOUBLE. 
  181.  *    
  182.  *
  183.  * Constants:
  184.  * The hexadecimal values are the intended ones for the following constants.
  185.  * The decimal values may be used, provided that the compiler will convert
  186.  * from decimal to binary accurately enough to produce the hexadecimal values
  187.  * shown.
  188.  *
  189.  */
  190.  
  191. #if defined(vax)||defined(tahoe)
  192. /*C0     =  4.1666666666666504759E-2    , Hex  2^ -4   *  .AAAAAAAAAAA9F0 */
  193. /*C1     = -1.3888888888865302059E-3    , Hex  2^ -9   * -.B60B60B60A0CCA */
  194. /*C2     =  2.4801587285601038265E-5    , Hex  2^-15   *  .D00D00CDCD098F */
  195. /*C3     = -2.7557313470902390219E-7    , Hex  2^-21   * -.93F27BB593E805 */
  196. /*C4     =  2.0875623401082232009E-9    , Hex  2^-28   *  .8F74C8FA1E3FF0 */
  197. /*C5     = -1.1355178117642986178E-11   ; Hex  2^-36   * -.C7C32D0A5C5A63 */
  198. static long C0x[]    = { _0x(aaaa,3e2a), _0x(a9f0,aaaa)};
  199. static long C1x[]    = { _0x(0b60,bbb6), _0x(0cca,b60a)};
  200. static long C2x[]    = { _0x(0d00,38d0), _0x(098f,cdcd)};
  201. static long C3x[]    = { _0x(f27b,b593), _0x(e805,b593)};
  202. static long C4x[]    = { _0x(74c8,320f), _0x(3ff0,fa1e)};
  203. static long C5x[]    = { _0x(c32d,ae47), _0x(5a63,0a5c)};
  204. #define C0    (*(double*)C0x)
  205. #define C1    (*(double*)C1x)
  206. #define C2    (*(double*)C2x)
  207. #define C3    (*(double*)C3x)
  208. #define C4    (*(double*)C4x)
  209. #define C5    (*(double*)C5x)
  210. #else    /* defined(vax)||defined(tahoe) */
  211. static double
  212. C0     =  4.1666666666666504759E-2    , /*Hex  2^ -5   *  1.555555555553E */
  213. C1     = -1.3888888888865301516E-3    , /*Hex  2^-10   * -1.6C16C16C14199 */
  214. C2     =  2.4801587269650015769E-5    , /*Hex  2^-16   *  1.A01A01971CAEB */
  215. C3     = -2.7557304623183959811E-7    , /*Hex  2^-22   * -1.27E4F1314AD1A */
  216. C4     =  2.0873958177697780076E-9    , /*Hex  2^-29   *  1.1EE3B60DDDC8C */
  217. C5     = -1.1250289076471311557E-11   ; /*Hex  2^-37   * -1.8BD5986B2A52E */
  218. #endif    /* defined(vax)||defined(tahoe) */
  219.  
  220. #define cos__C(z)    (z*z*(C0+z*(C1+z*(C2+z*(C3+z*(C4+z*C5))))))
  221.  
  222. extern int finite();
  223. extern double copysign(),drem();
  224.